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We derive first- and second-order piezoelectric coefficients for the zinc-blende III-V semiconduc¬ 
tors, {Al,Ga,In}-{N,P,As,Sb}. The results are obtained within the Heyd-Scuseria-Ernzerhof hybrid- 
functional approach in the framework of density functional theory and the Berry-phase theory of 
electric polarization. To achieve a meaningful interpretation of the results, we build an intuitive 
phenomenological model based on the description of internal strain and the dynamics of the elec¬ 
tronic charge centers. We discuss in detail first- and second-order internal strain effects, together 
with strain-induced changes in ionicity. This analysis reveals that the relatively large importance in 
the III-Vs of non-linear piezoelectric effects compared to the linear ones arises because of a delicate 
balance between the ionic polarization contribution due to internal strain relaxation effects, and the 
contribution due to the electronic charge redistribution induced by macroscopic and internal strain. 


I. INTRODUCTION 


Non-linear piezoelectricity and its importance in III-V 
semiconductors and nanostructures has been a topic of 
intense debate over the last few years.^^^The existence of 
non-negligible second-order piezoelectricity in GaAs and 
InAs was first proposed by Bester et al^ in the context of 
linear-response density functional theory (DET). It was 
soon revisited semi-empirically by Migliorato et alW us¬ 
ing Harrison’s bond orbital model and a strain-dependent 
Kleinman parameter ( calculated from DET. More re¬ 
cently, following Bester’s methodology, Beya-Wakata et 
alW have extended the calculation of first- and second- 
order coefficients to other zinc-blende (ZB) III-V bina¬ 
ries. Although failing for lower band gap materials due 
to the intrinsic limitations of the local density approxi¬ 
mation (EDA) , 1 ^ the approach by Bester et a/.® provides, 
for the higher energy gap III-V materials, ab initio pa¬ 
rameters without need to fit to experiment. Given that 
measuring piezoelectric (PZ) constants accurately is com¬ 
plicated, in particular when possibly large second-order 
effects are superimposed to the linear ones, the impor¬ 
tance of theoretical calculations cannot be overstated. It 
was indeed found by Beya-Wakata et al.^ and is con¬ 
firmed in the present work with the use of an improved 
energy functional, that second-order effects are large for 
the family of ZB III-V materials, in some cases dominat¬ 
ing over first-order piezoelectricity even for small strains 
(e.g. AlP, AlAs and InP). However, it remains unclear 
and hidden beneath the numbers what is the physical 
origin of this large second-order PZ effect, and what is 
the reason for the apparent lack of trends in the evolu¬ 
tion of the coefficients for the different compoundsWe 
try to address this question here by looking separately 
at the different factors which contribute to the PZ coef¬ 
ficients, investigating separately the linear and nonlinear 


contributions to the Kleinman parameter and to the ma¬ 
terial electronic response, as well as noting the influence 
of volume effects on the PZ coefficient values. 

The form of the second-order PZ polarization tensor for 
the ZB lattice, along with all the remaining PZ crystal 
classes, has been established by Grimmer based on sym¬ 
metry arguments.!^ The full expression for the PZ vector 
of a ZB lattice, including first- and second-order effects, 
reads P 


( e4 \ / eie4 \ 

65 I + Bii 4 I 6265 1 

ee / \ ^366 / 

( (^2 + £3)64 \ / 6566 \ 
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where 614 is the linear PZ coefficient, while Hi 14, H124 
and Hi 56 are the second-order PZ coefficients. The strain 
components are denoted by ej. We follow standard Voigt 
notation for all the symbols. 

Previous analysis by Beya-Wakata et al.^ found that 
Hii 4 and H124 show no obvious trend with compound. 
We note that the two terms in the first (x-) component 
of the polarization vector [Eq. Q] involving these two 
second-order PZ coefficients are proportional to 61 and 
^2 + ^3, respectively. However, it is often more useful and 
intuitive to describe the response of a material in terms 
of the hydrostatic strain Tr(6) = 61 + 62+63 and the 
biaxial strain 6b,i = ( 36 ^ — Tr( 6 ))/ 2 . Equation 0 can 
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then be rewritten 


P =ei 4 


eb,ie4 

■ ^2 I ^b,2^5 



( 2 ) 


with Ai and A 2 related to the more usual coefficients as 
follows: 


^1 = 


5ii4 + 2^124 


Ao 


2(5i14 — ^ 124 ) 


(3) 


While trends in ^114 and ^124 alone are difficult to ex¬ 
tract P we find that Ai and A 2 tend to decrease and to 
increase, respectively, with increasing unit cell volume. 
These opposing trends then account for the difficulty to 
identify any clear trends in Bn^ and B 124 ,. Overall, we 
show that, once the non-linear dependence on strain of 
the bond polarity, in terms of charge redistribution and 
internal strain effects, is taken into account, a complete 
explanation of the PZ effect, including its non linearity, 
can be provided from first principles alone. We present 
a consistent picture of ionicity in ZB materials, drawn 
from the well-established relation between the dynamics 
of Wannier function centers and electric polarization.^^ 
Although data on PZ coefficients is somewhat lacking 
for ZB semiconductors, more reliable measurements are 
available for the urtzite III-Ns. A very recent study of 
elastic and PZ properties of WZ GaN by Witczak et a/.^ 
serves as validation of the methodology employed here, 
which we already followed in Refs. [ 11 ] and [121 respec¬ 
tively. 

Recent work by Tse et relying on the Harrison 
model computed the effective PZ coefficient of GaAs and 
InAs using a third-order expansion in strain, which corre¬ 
sponds for a fourth-order expansion for the polarization. 
However, the authors neglected the quadratic contribu¬ 
tions arising from shear strain (related to the Bi^q coef¬ 
ficient) which, as we will show, are of the same order of 
magnitude as the other second-order contributions. 

The paper is organized as follows. We introduce our 
computational method in Sec. [Hj The main results are 
then presented in Sec. m This is followed by an anal¬ 
ysis of the results in Sec. m where we provide some 
general considerations in Sec. IV A[ before discussing the 
first-order and second-order coefficients in Secs. IIV 51 and 
IV C[ respectively. Finally, we summarize our conclu¬ 
sions in Sec.[Vl Some additional details of the calculation 
method used are presented in the Appendix. 


II. COMPUTATIONAL METHOD 

Equations 0 and 0 portray all the symmetries ap¬ 
plicable to piez oelectricity which are characteristic of the 
ZB lattice,^^^ and therefore form the basis of the target 


quantities being calculated here. The DFT calculations 
for geometry optimization and self-consistent wave func¬ 
tions were carried out using the H eyd-S cuseria-Ernzerhof 
(HSE) hybrid functional approaclP^ ^ wi thin the projec¬ 
tor augmented-wave (PAW) met ho as implemented 

in the VASP package.^I^ The screening parameter fi was 
fixed to 0.2 in our calculations and the mixing parame¬ 
ter a to 0.25 (these settings correspond to VASP’s HSE06 
version of the HSE functional). The plane-wave cutoff 
energy was set to 600 eV. The semicore d states of Ga 
and In were treated as valence states. A F-centered con¬ 
figuration was used for the /c-point mesh, with a variable 
number of divisions. The convergence of the PZ coeffi¬ 
cients was found to be slow with respect to the number 
of k points, as described further in the Appendix. Ta¬ 
ble |l] shows the calculated energy gaps and lattice pa¬ 
rameters for all the compounds considered. The calcu¬ 
lations of electric polarization were performed employ¬ 
ing the Berry-phase technique wit hin the context of the 
modern theory of polarization,I^^^^H^ as implemented by 
Martijn Marsman in VASP. We have successfully used a 
similar approach previously to study internal str ain an d 
electric polarization of wurtzite group-HI nitrides. 

In order to extract the PZ coefficients from the po¬ 
larization results, we have fit to second-order polynomi¬ 
als. We have performed calculations of for 5 different 
strain branches e^'^\ which read: 

= ( 0 , 0 , 0 ,/3,/3,^) ^ 

e(2) = (a,0,0,/3,0,0) ^ 

= (0,a,0,/3,0,0) ^ 
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Here the relation between unstrained and strained struc¬ 
tures is given by the strain transformation matrix, which 
for an arbitrary lattice vector at equilibrium = 
{Ra,l,Ra, 2 ,Ra, 3 ) IS givCU by 



i + ei f f ) 


^ Ra.,\ ^ 

K = 

f 1 + ^2 f 


Ra ,2 


\ 2 2 ^ ^ / V R<^i^ / 

We increment the a and P parameters in Eq. Q so that 
the Voigt components of the strain tensor follow the re¬ 
lations above and allow to calculate the corresponding 
PZ coefficients. Eor branch 9 data points per com¬ 
pound and per k mesh were calculated, up to \f3\ = 0.04. 
Eor the other branches (which involve two-dimensional 
fitting) 25 data points were used, with \a\ and \P\ up to 
0.02 and 0.04, respectively. We have verified that these 
strain ranges allow a good quality fitting for both linear 
and quadratic terms, as shown in the Appendix. Obvi¬ 
ously this leads to plenty of redundant information which 
is however useful in order to check consistency between 





















3 


TABLE L Direct band gap at T, indirect gap (when this is smaller than the direct one) and lattice parameter, for all the binary 
ZB III-V semiconductors, calculated from the HSE approach as explained throughout the text. A E-centered 10x10x10 
/c-point sampling is used for all of these calculations. 



Direct 

Al 

gap at r 

Ga 

(eV) 

In 


Al 

Indirect gap (eV) 

Ga 

In 


Lattice parameter (A) 

Al Ga In 

N 

5.61 

3.10 

0.56 

N 

4.58 

n/a 

n/a 

N 

4.365 

4.493 

4.991 

P 

4.21 

2.81 

1.43 

P 

2.31 

2.33 

n/a 

P 

5.471 

5.460 

5.904 

As 

2.86 

1.33 

0.41 

As 

2.10 

n/a 

n/a 

As 

5.687 

5.686 

6.116 

Sb 

2.23 

0.75 

0.31 

Sb 

1.76 

n/a 

n/a 

Sb 

6.188 

6.152 

6.563 


different calculations. These consistency checks are im¬ 
portant given the convergence issues with respect to the 
number of k points highlighted in the Appendix. They 
also allow to monitor the effect of symmetry reduction 
which we have shown to significantly affect the results of 
calculated elastic and structural properties of semicon- 
ductors.l^ 


III. RESULTS 

All our results for first- and second-order PZ coef¬ 
ficients of the 12 ZB binary compounds {Al,Ga,In}- 
{N,P,As,Sb} are summarized in Table where we 
present separately the ionic and electronic contributions 
to ei4, 5 ii4, 5 i 24 , ^156, and A2. A number of related 
results are presented in Fig. and will be discussed in 
more detail here and in the next section. While trends in 
Bii4 and ^124 are difficult to extract,!^ the magnitude of 
Ai shows a clear decrease with increasing unit cell vol¬ 
ume and group V row, as shown graphically in Fig. 

By contrast, A2 increases for fixed cation with increasing 
group V row. As mentioned above. Table |n| presents the 
ionic and electronic contributions to the different coef¬ 
ficients separately. It is readily clear from the numbers 
presented, in particular for 614, that the PZ response 
in ZB arises as the result of two much larger delicately 
balanced contributions: the ionic and electronic parts. 
These results are further discussed in Sec. HYB 

Second-order PZ coefficients for the ZB III-V family, 
except the nitrides, have already been obtained in the 
DFT frame by Beya-Wakata et al using the LDA.l^ Our 
results in Table [H] are overall very similar for the lin¬ 
ear PZ coefficients of the different compounds with wider 
band gap, but less so for the second-order coefficients, es¬ 
pecially for the narrow gap materials. This disagreement 
is not surprising since the LDA predicts a metallic state 
for some of these compounds (e.g. GaSb, InAs and InSb) 
which is not compatible with the existence of a macro¬ 
scopic electric polarization in the material. Two initial 
conclusions can be drawn from this comparison and from 
earlier work on WZ nitridesJ^^i) LDA and/or other DFT 
approximations can be sufficient to study linear PZ polar¬ 
ization when a positive band gap is correctly predicted 
and ii) when both schemes predict a positive gap, the 


main differences between standard DFT implementations 
and a hybrid functional approach arise from the disagree¬ 
ment in the predicted internal strain. This second con¬ 
sideration also accounts for the disagreement between 
LDA-DFT and HSE-DFT in the case of the spontaneous 
polarization of WZ nitrides, which arises directly from 
the disagreement regarding the WZ internal parameter 
itPXhe HSE approach has been shown to improve upon 
Kohn-Sham DFT for a wide range of semiconductors re¬ 
garding the description of not only ba nd gaps, b ut also 
lattice parameters and elastic properties ) ^^ * 1 5 | I6 J24] xhere- 
fore we expect also the description of internal strain to 
be superior when using the HSE functional. 


IV. DISCUSSION AND INTERPRETATION OF 
RESULTS 


In the following we discuss the results obtained in this 
work and interpret them using a simple physical ap¬ 
proach. Section [Tv A| briefly comments on separating the 
system into ionic and electronic parts and the significance 
of defining a reference starting point. Sections |IV B| and 
II V Cl then consider the results obtained for the first- and 
second-order piezoelectricity, respectively, and their in¬ 
terpretation within the context of the present model. 


A. General considerations 

One of the most striking features of piezoelectricity in 
the ZB HI-Vs is that, except for the highly ionic nitrides, 
the second-order PZ response is very large if compared to 
the first-order contribution. This assertion becomes ex¬ 
treme in the case of AlP and InP, for which the sign of the 
effective PZ coefficient as a function of hydrostatic strain 
^14 = ei4 -h AiTr(e) would reverse its sign for as little as 
0.2% tensile hydrostatic strain. As a matter of fact, 
Beya-Wakata et found a negative linear PZ coefficient 
for AlP, while our value in Table [TT| is positive. To picture 
the meaning of this change in sign of the effective PZ co¬ 
efficient one can imagine a ball-and-stick model in which 
cations and anions are represented by point charges of 
opposite sign. A positive 614 corresponds to a positively 
charged cation and a negatively charged anion, which is 
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TABLE IL All the different first- and second-order PZ coefficients of the ZB III-Vs, split into different contributions. See text 
for details. 


ei 4 (C/m^) 



Ionic part: (C/w?) 

Al Ga In 


Electronic part: 614 ® 
Al Ga 

(C/m^) 

In 


Total: ei 4 (G/m^) 

Al Ga In 

N 

-2.266 

-2.256 

-2.409 

N 

2.814 

2.622 

3.002 

N 

0.548 

0.366 

0.593 

P 

-1.553 

-1.437 

-1.503 

P 

1.568 

1.316 

1.518 

P 

0.014 

- 0.121 

0.016 

As 

-1.427 

-1.313 

-1.371 

As 

1.372 

1.108 

1.259 

As 

-0.055 

-0.205 

- 0.111 

Sb 

-1.236 

-1.169 

-1.183 

Sb 

1.141 

0.953 

1.021 

Sb 

-0.094 

-0.216 

-0.161 


Bii4 (C/m^) 



Ionic part 

Al 

. Dion 

’• 1^114 

Ga 

(C/m^) 

In 


Electronic part: Bfi^ 
Al Ga 

(C/m^) 

In 


Total: Bii 4 (G/m 
Al Ga 

In 

N 

13.08 

11.89 

11.40 

N 

-19.89 

-17.27 

-17.37 

N 

-6.81 

00 

CO 

id 

1 

-5.96 

P 

7.06 

6.36 

6.05 

P 

-9.08 

-7.59 

-7.59 

P 

- 2.02 

-1.23 

-1.54 

As 

6.38 

5.84 

5.43 

As 

-7.99 

-6.30 

-6.60 

As 

-1.61 

-0.99 

-1.17 

Sb 

5.25 

4.88 

4.37 

Sb 

- 6.01 

-5.18 

-5.00 

Sb 

-0.76 

-0.31 

-0.62 


Bi24 (C/m^) 



Ionic part: By^^ 

(C/m^) 


Electronic part: ^124 

(C/m^) 


Total: Bi 24 (G/m 



Al 

Ga 

In 


Al 

Ga 

In 


Al 

Ga 

In 

N 

9.78 

11.97 

9.74 

N 

-14.83 

-18.70 

-16.06 

N 

-5.04 

-6.73 

-6.32 

P 

5.98 

6.29 

5.99 

P 

-8.75 

-9.56 

-9.61 

P 

-2.76 

-3.27 

-3.62 

As 

5.41 

5.51 

5.30 

As 

- 8.00 

-8.72 

-9.62 

As 

-2.59 

-3.21 

-4.31 

Sb 

4.38 

4.29 

4.31 

Sb 

-6.37 

-7.06 

-8.36 

Sb 

-1.99 

-2.77 

-4.04 


Bi 56 (C/m^) 



Ionic part: B^i^q 

(C/m^) 


Electronic part: Bf^Q 

(C/m^) 


Total: Bi 56 (C/m^) 

1 


Al 

Ga 

In 


Al 

Ga 

In 


Al 

Ga 

In 

N 

6.32 

5.35 

3.56 

N 

-10.46 

-8.52 

-5.57 

N 

-4.15 

-3.18 

- 2.00 

P 

2.98 

3.13 

2.39 

P 

-4.42 

-4.51 

-3.41 

P 

-1.43 

-1.38 

- 1.02 

As 

2.82 

3.02 

2.27 

As 

-4.14 

-4.30 

-2.74 

As 

-1.32 

-1.28 

-0.46 

Sb 

2.15 

2.41 

2.02 

Sb 

-2.97 

-3.11 

-2.19 

Sb 

-0.82 

-0.70 

-0.16 


Ai (C/m^) 



Ionic part: 

(C/m^) 


Electronic part: Af^ 

(C/m^) 


Total: Al (G/m^ 

) 


Al 

Ga 

In 


Al 

Ga 

In 


Al 

Ga 

In 

N 

10.88 

11.95 

10.30 

N 

-16.52 

-18.22 

-16.50 

N 

-5.63 

-6.28 

- 6.20 

P 

6.34 

6.31 

6.01 

P 

- 8.86 

-8.90 

-8.94 

P 

-2.52 

-2.59 

-2.93 

As 

5.73 

5.62 

5.34 

As 

-7.99 

-7.91 

-8.62 

As 

-2.26 

-2.47 

-3.27 

Sb 

4.67 

4.48 

4.33 

Sb 

-6.25 

-6.44 

-7.24 

Sb 

-1.58 

-1.95 

-2.90 


A 2 (C/m^) 



Ionic part 

Al 

Ga 

(C/m^) 

In 


Electronic part: Af^ 
Al Ga 

(C/m^) 

In 


Total: A 2 (G/m^' 
Al Ga 

) 

In 

N 

2.20 

-0.05 

1.10 

N 

-3.38 

0.96 

-0.87 

N 

-1.18 

0.90 

0.24 

P 

0.72 

0.05 

0.04 

P 

- 0.22 

1.32 

1.34 

P 

0.49 

1.36 

1.38 

As 

0.64 

0.22 

0.09 

As 

0.01 

1.61 

2.02 

As 

0.65 

1.48 

2.09 

Sb 

0.58 

0.40 

0.04 

Sb 

0.24 

1.25 

2.24 

Sb 

0.82 

1.64 

2.27 
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FIG. 1. (Color online) Bar chart with the comparison of different internal strain parameters and PZ coefficients for the different 
III-V compounds. 


found to be the case for the nitrides, implying in the ball- 
and-stick model a valence electronic charge (to be added 
to the bare positive core charge) of < 3 electrons on the 
cation and > 5 electrons on the anion. 

This first-order result is then consistent with Harri¬ 
son’s model. In this model the bond polarity depends 
on the relative magnitude of the covalent interaction U 
and the difference C in the atomic orbital self-energies 
on the neighboring atomic sites. The total number of 
electrons on the cation site increases as the magnitude of 
the ratio oilljC increases. The large electronegativity 
of N atoms compared to the other group V elements then 
leads to a smaller ratio oiUjC for the nitrides, consis¬ 
tent with there being > 5 electrons on the N sites in the 
ball-and-stick model.^ We note however that the calcu¬ 
lated negative values of Ai are contrary to what would 
be predicted using the Harrison model, where the magni¬ 
tude of U decreases with increasing bond length, leading 
to a predicted flow of charge from the cation to the anion 
with increasing volume. However, the negative value of 
Ai describes the opposite behavior - transfer of charge 
from the anion to the cation with increasing volume. 

This disagreement in the direction of charge transfer 
may reflect factors omitted from the Harrison model, in¬ 
cluding electron correlation effects. Also, in the context 
of the modern theory of polarization, a representation 
in terms of point dipoles lacks rigorous support. Any 


attempt to achieve a discrete representation of the sys¬ 
tem as an ensemble of point charges requires to decouple 
ionic and electronic contributions, with the point elec¬ 
trons given by the Wannier centers of charge.^ For spin- 
polarized calculations the latter are assigned charge — 1, 
while in the case of spin degeneracy they are assigned 
charge —2. For the HI-Vs, the cation has core charge 
+3 and the anion has core charge +5. When semicore 
d states are considered for Ga and In, their core charges 
become +13; however in our discussion we will assume 
for simplicity that the corresponding extra valence elec¬ 
trons are localized around the core giving an effective core 
charge of +3. Then for ZB HI-Vs, in the spin-degenerate 
case, there are four centers of charge with charge —2 each, 
corresponding to the overall 8 valence electrons. This is 
depicted in Fig. for both the case where the cation 
is chosen to be located at the origin of the unit cell, in 
which we shall work preferentially since it is the stan¬ 
dard representation for ZB, and the case where the an¬ 
ion is located at the origin. Both representations must 
be (and are) equivalent for computing the polarization. 
For convenience, the four Wannier centers of charge can 
be further combined into their center of mass, which is 
then assigned charge —8. 
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Cation-centered origin 



Anion-centered origin 



FIG. 2. (Color online) Zinc-blende lattice primitive unit cells 
showing the underlying cubic symmetry in both cation-at- 
origin (top) and anion-at-origin (bottom) conventions. Large 
(red) balls represent cations, with core charge -h3, while small 
(blue) balls represent anions, with core charge -h5. Solid balls 
are cores contained in the unit cell of the hgure and shaded 
balls are cores contained in neighboring unit cells. The dark 
dots placed along the bonds represent (spin-degenerate) Wan- 
nier centers of charge, with charge —2. 


B. First-order internal strain and piezoelectricity 


The ionic dipole moment in the cation-at-origin rep¬ 
resentation is given by the relative position of the anion 
with respect to the cation, which is kept fixed: 

d-- = 5e[(l + e)ro+t]. (6) 

Here, fq is the position of the anion before strain, and 
t is the anion’s internal strain, which for ZB is given, to 


first order^ in terms of the Kleinman parameter (^1111211111 


t(i) 




(7) 


The appearance of a strain-dependent part in Eq. 
is a consequence of the choice of origin to compute the 
dipole moment, and was already discussed in the seminal 
paper of Vanderbilt and King-Smith.^^ Alternatively, if 
we had placed the origin to compute the dipole moment 
at a position exactly 5/8 along the vector connecting the 
cation and the anion (keeping the same convention for the 
unit cell), strain dependence would have been removed: 


d^°^ =5e 


_ .3ro 3t 


T 3e 


(1 + ^) 


-5ro 


_15e 


( 8 ) 


reaching the result that the ionic dipole moment in ZB 
depends only on internal strain. Therefore any depen¬ 
dence on strain of this value must arise from the de¬ 
pendence of t on strain. Since the result must be in¬ 
dependent of the choice of origin, the ionic dipole mo¬ 
ment can be obtained as in Eq. ^ and then the 
strain-dependent artifact subtracted. Therefore, in the 
cation-at-origin convention, using Eqs. © and 0 we 
can express the linear ionic PZ coefficient as 


efl = -^C, (9) 

where the term appears after dividing by the unit cell’s 
volume V = a^/4. 

Conversely, within this same convention, the linear 
electronic PZ coefficient can be expressed as 

efi = (10) 


which in turn serves as the definition for the “elec¬ 
tronic Kleinman parameter” Note that the strain- 

dependent artifact for the electronic part arises in the 
same fashion as in Eq. 0 if Wannier centers are used, or 
arises as the so-called “quantum of polarization” if the 
computation is done directly from the Berry phase.^^ 
Einally the total linear PZ coefficient is given by 

ei4 = e-(r + e?!f = 4(8C‘^-5C). (11) 


Within the anion-at-origin convention, the ionic and elec¬ 
tronic contributions would have been given by 3e (/and 
8e(C®^® — C)/a^, respectively, thereby, as expected, leaving 
the total PZ coefficient unaltered. 

Equation 0 provides an extremely useful and intu¬ 
itive insight to the origin of the PZ effect in ZB III-V 
semiconductors. A net dipole moment in the unit cell 
appears as a consequence of the loss of centrosymmetry 
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and can be expressed in terms of the internal strain pa¬ 
rameters (see also Ref. [ 12 ] for an equivalent treatment 
in wurtzite). The magnitude of this PZ effect is then 
given by the ability (or inability) of the Wannier cen¬ 
ters of charge to follow the movement of the anion. This 
idea was already highlighted by Vanderbilt“t/z-e proper 
piezoelectric response is identically zero for a homoge¬ 
neous [strain] deformation of both the ionic positions and 
the Wannier centers^\ Rearranging Eq. 0 as follows 
allows to emphasize this idea even further: 


TABLE III. Kleinman parameter and ratio of the “electronic 
Kleinman parameter” with respect to ^ for all the III-Vs. 



Kleinman parameter ^ 


Ratio 

/c 


Al 

Ga 

In 


Al 

Ga 

In 

N 

0.539 

0.568 

0.749 

N 

0.776 

0.726 

0.779 

P 

0.580 

0.535 

0.654 

P 

0.631 

0.572 

0.631 

As 

0.576 

0.530 

0.640 

As 

0.601 

0.527 

0.574 

Sb 

0.591 

0.552 

0.636 

Sb 

0.577 

0.510 

0.539 


< foC 


ei 4 = ( S— -5 


Aele 


( 12 ) 


where the ratio characterizes the ability of the 

Wannier centers of charge to follow the anion core, 
^eie/^ _ 2 corresponds to these centers following the core 
exactly, that is, it depicts an ideal ionic picture where the 
charge of the anion can be effectively expressed as its ox¬ 
idation number (—3 in this case). = 0 corresponds 

to the Wannier centers fully ignoring the movement of 
the anion core and remaining at their original positions 
in the lattice (plus a strain transformation). This situa¬ 
tion would interestingly lead to a larger PZ constant, of 
negative sign, indicating reversed ionicity, i.e. the centers 
of charge follow the cation movement exactly. Eor III-Vs, 
charge balance is achieved when = 5/8 = 0.625, 

in which case the PZ response of the material is zero. 
The Kleinman parameter for the different compounds is 
shown in Table [in| (also in graphical form in Eig. to- 
gether with the ratio where is obtained from 

Eq. ( P!q| ). All the compounds with C®^®/C > 0.625 present 
regular ionicity and a positive PZ coefficient (basically 
the nitrides), while the rest present reversed ionicity, with 
a net positive charge on the anions. The proximity of 
the value of fo 0.625 for many of these materials 

then explains the importance of the quadratic terms 
computing the PZ response of the III-Vs, and even 
determining its sign. 

Eor comparison, we have computed the numbers for 
extreme case of ionicity, ZB NaE, a I-VII compound for 
which C = 0.921 and C®^®/C = 0.982, that is the Wannier 
centers of charge follow the anionic core almost exactly. 
A Kleinman parameter C, close to 1 is also indicative of 
high ionicity because it implies the response of the asym¬ 
metric unit to strain is to preserve the tetrahedral bond 
lengths rather than the bond angles 

We see from Table |III| and Eig. that the Kleinman 
parameter has a relatively small variation across the dif¬ 
ferent III-V compounds considered, although its value 
does tend to be slightly higher in the In-V compounds, 
reflecting that the bond-bending force constants are rel¬ 
atively weaker in the In-V compounds compared to the 
Al- and Ga-containing compounds. We also note that 
the PZ coefficient of NaE (and other I-VIIs) is not given 

by Eq. (12), but by 614 = -7) instead. This 


analysis alTows to establish the maximum linear PZ co¬ 
efficient that can exist for the different families of binary 


ZB compounds as a function of their lattice parameter as 


oinax _ Xe 

-14 


= where X is the absolute value of the anion’s 
oxidation number. In practice, 614 is only close to this 
value for the I-VII compounds. 

Another important consideration that becomes imme¬ 
diately clear from Eq. (12) is the volume effect present 


in piezoelectricity, whereby the lattice parameter enters 
the expression of the PZ coefficient as 1/a^. This factor 
varies for the ZB III-V compounds between 0.023 
for InSb and 0.052 for AIN, so that at equal dipole 
moment per atom pair, this volume effect would enhance 
the PZ coefficient of AIN by more than a factor of 2 com¬ 
pared to InSb. 


C. Second-order internal strain and piezoelectricity 

The second-order internal strain parameters in ZB, 
which we call Uij , can be obtained in a similar way to the 
Kleinman param eteJ^ and the wurtzite internal strain 
parameters,im^ expanding the expression of the inter¬ 
atomic energy to third-order in strain and evaluating the 
different terms: 


in 

=az/i4 

' 6164 ^ 


in 

^2^5 

+ az/24 



\ ^3^6 ) 



^ (^2 + ^3)^4 ^ 

(ei + 63)65 


az /56 



(13) 


Equation ( [l^ can also be expressed in terms of hydro¬ 
static and biaxial strain: 


=az/iTr(6) 


az /56 



az/2 



(14) 


It is unsurprising that Eqs. and ( p^ resemble the ex¬ 
pressions of the second-order PZ polarization [cf. Eqs. 0 
and respectively]. We will not report their values for 
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FIG. 3. (Color online) Evolution as a function of hydrostatic and biaxial strain of the strain-dependent Kleinman parameter 
Cl and Cf^/Ci ratio for all the ZB III-V binaries. The distance between adjacent contour lines is 0.125 in the color-coded scale. 


all the III-Vs here. They can however be computed from 
the second-order ionic PZ coefficients in Table nn in a 
similar manner to the linear one: 


Z/14 


4 5e 



Dion 

-^156 

5e 




4 5e 


(15) 


and 


_a2 

“T 5 ^ 


J^2 




(16) 


where the term involving in the different z/s comes 
from the change in volume due to the hydrostatic strain. 
LookiM at the coefficients for in 

Table and taking the scaling factor from Eq. (15) 
into account, it can be seen that the second-order internal 
strain shows only a relatively small variation between the 
different binaries, as was also the case for the Kleinman 
parameter (Table [nT| and the bar charts in Fig. [^. 

The equivalent expressions for the center of mass of 


the Wannier centers are: 


ele 

^14 


ele 

^56 


_a2 + efi 

■4 -8e ’ 

^2 Dele 
-^156 

"4 -8e ’ 


ele 

^24 


+ > 


jele 


4 -8e 


and 


= 


a? 


14 


-8e 


, ,ele _ 
Z /9 — 




(17) 


(18) 


Using the expression for the second-order internal strain, 
Eq. (13), together with the first order one, Eq. ([^, allows 
to construct strain-dependent analogues of C and as 
a function of hydrostatic and biaxial strain: 


Ci(6) = C - 4z/iTr(6) - 4z/2eb,i, 

^ele(^) ^ ^ele _ 


(19) 


with 


, ,ele 


= 

,ele 


Z^14 + 2z/24 


Z^2 = 


2z/i4 — 2z^24 


G4 


2vfi 


,ele 




( 20 ) 
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It can be seen from the values of and in Table [n| 
that the magnitude of Ci{€) tends to increase with de¬ 
creasing volume (Tr(e) < 0), as would be expected from 
the bonds becoming more rigid with decreasing volume. 
By contrast, Ci{^) has a much weaker dependence on bi¬ 
axial strain, eb,i. 

We can in addition combine the two terms in Eq. (19) 
to rewrite Eq. (12) in the following form: 


e 


i _ 

14 — 


eCi 



( 21 ) 


extremely delicate balance between these two large com¬ 
ponents. 

We have built a phenomenological model of the PZ 
effect in the III-Vs that allows a straightforward inter¬ 
pretation of first- and second-order effects in terms of 
changes in internal strain and ionicity. While we find 
changes in internal strain to be largely dominated by hy¬ 
drostatic strain alone, changes in ionicity are found to be 
also strongly sensitive to biaxial strain. These changes 
are more pronounced as the unit cell size increases: they 
are smallest for AIN and largest for InSb. 


The second-order contribution to the PZ response then 
depends both on the evolution of Q and of Cf^/Ci- h 
Cf^/Ci was independent of strain then it would be im¬ 
possible for 6^4 to change sign with increasing strain. In 
Fig-i we plot the variation of the strain-dependent pa¬ 
rameter Q (left) and the strain-dependent ratio Cf^/Q 
(right) for all the ZB III-Vs, as a function of both hydro¬ 
static strain (change in volume) and biaxial strain. We 
see in the left hand figure for Q that contours of con¬ 
stant Q are generally close to vertical, highlighting that 
the hydrostatic strain has a considerably stronger effect 
than biaxial strain on the effective Kleinman parameter. 
The right figure elucidates how strain induces a transition 
(denoted by the white region in the color maps, for which 
Cf^/Ci = 5/8) from a negatively charged to a positively 
charged anion in all the compounds except for the ni¬ 
trides (only GaN at very high tensile hydrostatic strain 
sees this transition for the range shown). Eor many of 
the compounds this transition occurs remarkably close to 
the zero strain region, which highlights the importance 
of the second-order PZ coefficients. We also note that 
the contours of constant Cf^/Q are close to vertical for 
the III-nitrides, becoming more tilted toward the bottom 
right hand corner of the figure, reflecting the trend pre¬ 
viously noted that Ai tends to decrease and A 2 tends to 
increase with increasing compound’s unit cell volume. It 
is worth noting also that the nitrides retain their ionic 
character, a fact that is also reflected in their relatively 
constant Born effective charge when compared to other 
III-V compounds.!^ 


V. SUMMARY 

In summary, we have calculated hybrid-functional fir st¬ 
and second-order PZ coefficients for the whole family of 
zinc-blende III-V semiconductors. The hybrid-functional 
approach allows a more accurate evaluation of structural 
parameters and, more particularly, of band gaps. Other 
DET approximations predict a negative band gap for 
some of these compounds which prevents the calculation 
of a meaningful electric polarization. We have reported 
the PZ coefficients split into their ionic and electronic 
contributions and have shown how the large second-order 
corrections to the first-order coefficients arise from the 
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Appendix: Calculation details 

The convergence of the electric polarization, and there¬ 
fore of the PZ coefficients, is rather slow with respect to 
the number of k points used in the calculation. This 
becomes problematic with the highly demanding HSE 
functional, where typical computational times are one 
to two orders of magnitude higher than standard imple¬ 
mentations of Kohn-Sham DET. The effect of the finite 
k sampling on the data affects both the slope at zero 
strain and the curvature of the polarization, that is first- 
as well as second-order PZ coefficients. Even for a large 
number of k points (10 x 10 x 10 ) the derived values are 
not fully converged. We have found, however, a clear 
asymptotic behavior as the number of k points grows, as 
can be inferred from Fig. I^for GaSb under homogeneous 
shear deformation (64 = 65 = ee) for different k grids. 
The non-zero value of the polarization at zero strain in 
this case can be observed to be a numerical error due 
to finite k sampling.l^ However, the error is consistent 
for the different data points within the same series and 
therefore the extracted PZ coefficients are independent of 
it. In this context, the values of the extrapolated PZ co¬ 
efficients have been obtained assigning different weights 
to the coefficients obtained with different precision, with 
the lowest weight assigned to those obtained with 4x4x4 
k meshes and the highest weights corresponding to the 
denser lOxlOxlO/c meshes. Intermediate 6 x 6 x 6 
and 8 x 8 x 8 meshes have also been considered. In 
the LDA, an additional number of k points is specifi¬ 
cally required for materials for which a (wrong) negative 
band gap is predicted, as is the case of GaSb, InN, InAs 
and InSb. Beya-Wakata et a/., for instance, used a very 
dense 12 x 12 x 12 /c mesh to obtain their LDA results.® 
In our HSE calculations all the HI-Vs studied retain a 
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FIG. 4. (Color online) Convergence of the electric polarization 
values in GaSb as a function of homogeneous shear strain, for 
different Brilloin zone samplings. 


nearly correct (as compared to experiment) positive gap, 


as shown in Table [T[ This renders the matter of /c-point 
density less critical. 

To compensate for finite /c-mesh sampling, we assume 
that the absolute value of the numerical error s in either 
eij or Bijj^ due to finite Brilloin zone integration decays 
exponentially with the number of k points Nk used in 
the calculation: 

\s\=ae~^^^^ a^l3>0. (^- 1 ) 

For the linear coefficient 614 , we assume the error is 50% 
for the 4x4x4 mesh and 1% for the 10 x 10 x 10 mesh. For 
the non-linear coefficients Bij^ these errors are assumed 
to be 100% and 5%, respectively. These numbers yield 
a = 0.653 and P = 4.18 X 10 ^ for the linear coefficients 
and a = 1.227 and p = 3.2 x 10“^ for the non-linear 
ones. Then the coefficients are obtained, for each strain 
branch, as a weighed average of the four available values 
(corresponding to k-meshes 4x4x4, 6 x 6 x 6 , 8 x 8 x 8 
and 10 X 10 X 10), with weights w{Nk) = l/s‘^. The 
final coefficient presented in this paper is obtained as 
the average of these values between the different strain 
branches used. 
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